#include <iostream>
#include <fstream>
#include <string>
#include <sstream>
#include <list>
#include <cmath>

#include "Meshes/StructuredMesh.hpp"
#include "Meshes/Uniform.hpp"
#include "Utils/inout.hpp" 

#include "ParTrack/ParticleTracking.hpp"
#include "ParTrack/Linear.hpp"
#include "ParTrack/RungeKutta4.hpp"

using namespace std;
using namespace Tuna;

typedef TunaArray<double, 3>::tiny floatTinyArray_t;
typedef TunaArray<double, 3>::huge ScalarField3D;
typedef TunaArray<double, 1>::huge ScalarField1D;

void writeToFile_DX(list<floatTinyArray_t> &, int, string);
void writeToFile_GNU(list<floatTinyArray_t> &, int, string);
void writeToFile_VTK(const ScalarField3D&, std::string, std::string);
void readFromFile_DX_3D(ScalarField3D&, ScalarField3D&, ScalarField3D&, 
			int, string, int w = 5); 

timer crono;

int main () 
{
  double length_x, length_y, length_z, Z0, dt, pX1, pY1, pX2, pY2, tol, dist;
    double s0, cD;
    int num_nodes_x, num_nodes_y, num_nodes_z;
    int NpX, NpY, frec, Nstart, Nend, jump, Nxf, Nyf, Nzf;

    double etime = 0, factNormal;
// -----> Reading data from the input file
    std::ifstream input_cin ("inputPoincare");
    input_cin >> length_x       // Length of the 3D domain
	      >> length_y       // Length of the 3D domain
	      >> length_z       // Length of the 3D domain
	      >> Z0             // Thickness
	      >> num_nodes_x    // Number of nodes in X 
	      >> num_nodes_y    // Number of nodes in Y
	      >> num_nodes_z    // Number of nodes in Z
	      >> pX1   // X-coordinates for the plane
	      >> pY1  // X-coordinates for the plane
 	      >> pX2   // X-coordinates for the plane
	      >> pY2  // X-coordinates for the plane
	      >> NpX
	      >> NpY
	      >> Nstart
	      >> Nend
	      >> jump
	      >> dt
	      >> s0
	      >> cD
	      >> tol
	      >> Nxf
	      >> Nyf
	      >> Nzf
	      >> frec
	      >> factNormal;
    
    input_cin.close();


// -----> Mesh for velocity
    StructuredMesh< Uniform<double, 3> > mesh(length_x, num_nodes_x, 
					      length_y, num_nodes_y,
					      length_z, num_nodes_z);
    floatTinyArray_t deltas = mesh.getDeltas();
    mesh.print();
    
// -----> List for points storage
    list<floatTinyArray_t> linec;
    list<floatTinyArray_t>::iterator iter;

    floatTinyArray_t aux; 
    floatTinyArray_t x_i_n;     // Coordinates of x(i) at step n
    floatTinyArray_t x_i_np1;   // Coordinates of x(i) at step n+1
    floatTinyArray_t x_im1_n;   // Coordinates of x(i-1) at step n
    floatTinyArray_t x_im1_np1; // Coordinates of x(i-1) at step n+1
    floatTinyArray_t p1, p2, p3, p4;
    double deltaPX = (pX2 - pX1) / (NpX - 1);
    double deltaPY = (pY2 - pY1) / (NpY - 1);
    double deltaL = deltaPX;
    double s0_deltaL = s0 * deltaL;
    double deltaL0 = deltaL;

    cout << "\n +----- Line construction -----+ ... ";

    list<double> s, tao, rho, delta;
    list<double>::iterator iS, iT, iR, iD;

    for(int j = 0; j < NpY; j++) {
      for(int i = 0; i < NpX; i++) {
	aux(0) = pX1 + i * deltaPX;
	aux(1) = 0.5;
	aux(2) = 0.5 + j * deltaPY;
	linec.push_back(aux);
      }
    }


    /***
    for(int i = 0; i < NpX - 1; i++) {
	s.push_back(s0);
	tao.push_back(0.0);
	rho.push_back(1.0);
	delta.push_back(deltaL0);	
    }
    /***/
    cout << " +----- DONE -----+\n";

    writeToFile_DX(linec, 0, "./DataMix3D_01L/plan.");

    cout << "\n +--------------------------+"
	 << "\n +      Begining track      |"
	 << "\n +--------------------------+\n\n";

// -----> Velocity field at time t
    ScalarField3D u(mesh.getExtentNodes());       
    ScalarField3D v(mesh.getExtentNodes());
    ScalarField3D w(mesh.getExtentNodes());

// -----> Velocity field at time t + dt
    ScalarField3D un(mesh.getExtentNodes());       
    ScalarField3D vn(mesh.getExtentNodes());
    ScalarField3D wn(mesh.getExtentNodes());

    ParticleTracking< RungeKutta4< Linear < double, 3> > > particle;


    /***
    //    ScalarField1D s(NpX-1), tao(NpX-1);
    ScalarField3D c(Nxf, Nyf, Nzf);
    double sigmaXx, sigmaXy, sigmaXz, sigmaZx, sigmaZy, sigmaZz;
    double coordx, coordy, coordz, midx, midy, midz, normalx, normaly, normalz;
    double aux1, aux2, aux3, normSigX, normSigZ;
    double dxf = length_x / (Nxf - 1);
    double dyf = length_y / (Nyf - 1);
    double dzf = length_z / (Nzf - 1);
    double area;
    ***/

    cout << " +----- Time series -----+\n"
	 << " | Nstar = " << Nstart << "\t Nend = " << Nend 
	 << "\t jump = " << jump 
	 << "\n +-----------------------+ ";
// -----> Main loop

    int pausa;

    /***/
    for(int step = Nstart; step < Nend; step += jump) {      

      cout << "\n\n ----> Step = " << step;

      cout << "\n   +----- Tracking ... ";
      crono.tic();
// -----> Read the velocity at time t from a file
      readFromFile_DX_3D(u, v, w, step, "../../06ConvDiff/DataMix3D_01/velc.");
//      InOut::readFromFile_1(u, v, step, "../Velocity/Data01/velc.");

// -----> Read the velocity at time t + dt from a file
      readFromFile_DX_3D(un, vn, wn, step + jump, "../../06ConvDiff/DataMix3D_01/velc.");
//      InOut::readFromFile_1(un, vn, step + jump, "../Velocity/Data01/velc.");

//      area = 0.0;
// -----> Track every particle on the line
      iter = linec.begin(); // Pointer to initial point in the line x(i-1)^n
      x_im1_n = *iter;      // Contents of initial point in the line x(i-1)^n

// Track first point to the next position x(i-1)^(n+1) 
      *iter = particle.track(*iter, deltas, u, v, w, un, vn, wn, dt);
      x_im1_np1 = *iter; 
      iter++;  // Advance iter to the second point in the line

      for( ; iter != linec.end(); ++iter) {
	x_i_n = *iter; // Current point  x(i)^n
	x_i_np1 = particle.track(*iter, deltas, u, v, w, un, vn, wn, dt); 

/****
/// Calculate the distance between consecutive tracked points 
/// x(i-1)^(n+1) and x(i)^(n+1)
	dist = sqrt( ( x_i_np1(0) - x_im1_np1(0) ) * 
		     ( x_i_np1(0) - x_im1_np1(0) ) +
		     ( x_i_np1(1) - x_im1_np1(1) ) * 
		     ( x_i_np1(1) - x_im1_np1(1) ) +
		     ( x_i_np1(2) - x_im1_np1(2) ) * 
		     ( x_i_np1(2) - x_im1_np1(2) ) );
/// if dist exceeds the tolerance, adds a point between points: 
/// x(i)^n and x(i-1)^n
	
	if ( dist > tol ) {
	  aux(0) = ( x_i_n(0) + x_im1_n(0) ) * 0.5;
	  aux(1) = ( x_i_n(1) + x_im1_n(1) ) * 0.5;
	  aux(2) = ( x_i_n(2) + x_im1_n(2) ) * 0.5;
	  linec.insert(iter,aux);
	  iter--; iter--;

	  //---> Estaba comentado	  
	  s.insert(iS, *iS); // Striation thickness will keep its value for the new two intervals.
	  iS--;iS--; // Rewind iS one position, to the previously inserted striation thickness

	  tao.insert(iT, *iT); // Dimensionless time will keep its value for the new two intervals.
	  iT--;iT--; // Rewind iT one position, to the previously inserted tao.

	  delta.insert(iD, dist * 0.5);
	  iD--; iD--;	  
	  //---> Estaba comentado	  

	} else {

	  x_im1_n   = x_i_n;
	  x_im1_np1 = x_i_np1;
	  *iter     = x_i_np1;
	  
	  //---> Estaba comentado	  
	  *iS =  s0 * deltaL0 / dist; // Striation thickness: s_i = s0 * delta_0 / delta_i
	  *iT = *iT + dt * cD / ( (*iS) * (*iS) ); // Calc. dimensionless time

	  deltaL0 += *iS * sqrt(1 + 4 * (*iT));
	  //---> Estaba comentado	  

	  area += dist * Z0; // Accumulate the length of the strip	 	  

	}
/****/
      }

      cout << " Elapsed time = " << crono.toc();

      if( !(step % frec) ) {
	cout << "\n   +----- Writing file ... ";
	crono.tic();
	//	writeToFile_DX(linec, step, "./DataMix3D_01L/plan.");
	writeToFile_DX(linec, step, "/tmp/DATOS/DataMix3D_01/plan.");
	cout << " Elapsed time = " << crono.toc();
      }
      cout << "\n   +----- Diffusion calculation ... " << flush;
    
      /* *
// ------> 1) Calc. striation thickness and dimensionless time
      crono.tic();

      iS = s.begin(); iT = tao.begin();       
      iter = linec.begin();

      //deltaL = 0.0;
      for( ; iS != s.end(); ++iS, ++iT) {
	p1 = *iter; iter++; p2 = *iter; iter++; p3 = *iter; 
	*iS = s0 * length_z * deltaL0 / ( sqrt( (p1(0) - p3(0)) * (p1(0) - p3(0)) +
						(p1(1) - p3(1)) * (p1(1) - p3(1)) +
						(p1(2) - p3(2)) * (p1(2) - p3(2)) ) * 
					  sqrt( (p1(0) - p2(0)) * (p1(0) - p2(0)) +
						(p1(1) - p2(1)) * (p1(1) - p2(1)) +
						(p1(2) - p2(2)) * (p1(2) - p2(2)) ) );
	*iT = *iT + dt * cD / ( (*iS) * (*iS) );

	//deltaL += *iS * sqrt(1 + 4 * (*iT)); 
	
      }    

      //deltaL /= (NpX - 1);

      /* *     
// ------> 2) 
      c = 0;

      for(int k = 0; k < Nzf; ++k)
	for(int i = 0; i < Nxf; ++i)
	  for(int j = 0; j < Nyf; ++j) {
	    iS = s.begin(); iT = tao.begin(); iter = linec.begin();
	    for( ; iS != s.end(); ++iS, ++iT) {	  
	      
	      p1 = *iter; iter++;
	      p2 = *iter; iter++;
	      p3 = *iter;

	      sigmaXx = p3(0) - p1(0); 
	      sigmaXy = p3(1) - p1(1); 
	      sigmaXz = p3(2) - p1(2); 
	      normSigX = sqrt ( sigmaXx * sigmaXx + 
				sigmaXy * sigmaXy + 
				sigmaXz * sigmaXz );
	      sigmaXx /= normSigX; sigmaXy /= normSigX; sigmaXz /= normSigX;

	      sigmaZx = p2(0) - p1(0); 
	      sigmaZy = p2(1) - p1(1); 
	      sigmaZz = p2(2) - p1(2);

	      normSigZ = sqrt ( sigmaZx * sigmaZx + 
				sigmaZy * sigmaZy + 
				sigmaZz * sigmaZz );
	      sigmaZx /= normSigZ; sigmaZy /= normSigZ; sigmaZz /= normSigZ;

	      normalx = sigmaZy * sigmaXz - sigmaZz * sigmaXy;
	      normaly = sigmaZz * sigmaXx - sigmaZx * sigmaXz;
	      normalz = sigmaZx * sigmaXy - sigmaZy * sigmaXx;

	      normSigZ = sqrt ( normalx * normalx +
				normaly * normaly +
				normalz * normalz );
	      normalx /= normSigZ; normaly /= normSigZ;  normalz /= normSigZ;

	      coordx = dxf * i;
	      coordy = dyf * j;
	      coordz = dzf * k;

	      midx = p3(0); 0.5 * ( p2(0) + p3(0) ); 
	      midy = p3(0); 0.5 * ( p2(1) + p3(1) );	  
	      midz = p3(0); 0.5 * ( p2(2) + p3(2) );	  

	      aux1 = (coordx - midx) * sigmaXx + 
		     (coordy - midy) * sigmaXy + 
		     (coordz - midz) * sigmaXz;

	      aux2 = (coordx - midx) * sigmaZx + 
		     (coordy - midy) * sigmaZy + 
		     (coordz - midz) * sigmaZz;

	      aux3 = (coordx - midx) * normalx + 
		     (coordy - midy) * normaly + 
		     (coordz - midz) * normalz;

	      aux1 = aux1 * aux1 / (deltaL * deltaL);
	      aux2 = aux2 * aux2 / (deltaL * deltaL);
	      aux3 = aux3 * aux3 * factNormal/ ( (*iS) * (*iS) * (1.0 + 4.0 * (*iT)) );

	      c(i,j,k) += exp(-aux1-aux2-aux3) / (1.7726 * sqrt( 1.0 + 4.0 * (*iT)));
	      //	      cout << exp(-aux1-aux3) << " ";
	    }
	  }

      cout << " Elapsed time = " << crono.toc();
      
      cout << "\n   +----- Writing file ....";
      crono.tic();

      
      c /= max(c);	
      InOut::writeToFile_DX(c, step, "./DataMix3D_01L/conc.", dxf, dyf, dzf); 	    
      cout << " Elapsed time = " << crono.toc();
      /* */
    }

    /* */
    
    cout << "\n End ... \n\n";

    return 0;
}

void writeToFile_DX(list<floatTinyArray_t> &linec, int n, string filename)
{

    int Npoints = linec.size();
    list<floatTinyArray_t>::iterator iter;

    ostringstream extnumber;
    extnumber.width(5);
    extnumber.fill('0');
    extnumber << n;
    string snumber = extnumber.str();
    filename += snumber;
    ofstream file (filename.c_str());  

    //    cout << "writing " << filename << endl;
    
    file << " object 1 class array type float rank 1 shape 3 items "
	 << Npoints << " data follows\n";
    for(iter = linec.begin(); iter != linec.end(); ++iter) {
      file << (*iter)(0) << "   " << (*iter)(1) << "   " << (*iter)(2) << "   ";
    }
    file << endl;

    file << "object 2 class array type int rank 1 shape 2 items "
	 << Npoints - 1 << " data follows\n";
    for(int i = 0; i < Npoints - 1; ++i)
      file << i << "   " << i + 1 << "   ";
    file << endl;

    file << "attribute \"element type\" string \"lines\"\n"
	 << "attribute \"ref\" string \"positions\"\n";
    file << "object 3 class array type float rank 0 items "
	 << Npoints << " data follows\n";
    for(int i = 0; i < Npoints; ++i)
      file << 1 << "\t";
    file << endl;
 
    file << "attribute \"dep\" string \"positions\"\n";
    file << "object \"irregular positions regular connections\" class field\n"
	 << "component \"positions\" value 1\n"
	 << "component \"connections\" value 2\n" 
	 << "component \"data\" value 3\n" 
	 << "end\n" ;
    file.close();

}

void writeToFile_GNU(list<floatTinyArray_t> &linec, int n, string filename)
{

    int Npoints = linec.size();
    list<floatTinyArray_t>::iterator iter;

    ostringstream extnumber;
    extnumber.width(5);
    extnumber.fill('0');
    extnumber << n;
    string snumber = extnumber.str();
    filename += snumber;
    ofstream file (filename.c_str());  

    cout << "writing " << filename << endl;
    
    for(iter = linec.begin(); iter != linec.end(); ++iter)
	file << (*iter)(0) << "\t" << (*iter)(1) << "\n";
    file << endl;

    file.close();

}

void writeToFile_VTK(const ScalarField3D& data, 
		     std::string name, std::string variable)
{
    int Nvtk_x = data.rows();
    int Nvtk_y = data.cols();
    int Nvtk_z = 1;
    double dx = 1. / (Nvtk_x - 1);
    double dy = 1. / (Nvtk_y - 1);
    double dz = 1.;
    int Ntotal = Nvtk_x * Nvtk_y * Nvtk_z;

    std::ofstream file (name.c_str());        
    file << "# vtk DataFile Version 2.0 \n"
	 << "Lid_driven streamfuntion-vorticity \n"
	 << "ASCII \n"
	 << "DATASET RECTILINEAR_GRID \n" 
	 << "DIMENSIONS " 
	 << Nvtk_x << " " << Nvtk_y << " " << Nvtk_z;

    file << "\nX_COORDINATES " << Nvtk_x << " float \n";
    for(int i = 0; i < Nvtk_x; ++i) {
	file << i * dx << " ";
	if ( !(i % 6) ) file << "\n";
    }

    file << "\nY_COORDINATES " << Nvtk_y << " float \n"; 
    for(int i = 0; i < Nvtk_y; ++i) {
	file << i * dy << " ";
	if ( !(i % 6) ) file << "\n";
    }

    file << "\nZ_COORDINATES " << Nvtk_z << " float \n"; 
    for(int i = 0; i < Nvtk_z; ++i) {
	file << i * dz << " ";
	if ( !(i % 6) ) file << "\n";
    }
	
    file << "\nPOINT_DATA " << Ntotal
	 << "\nSCALARS " << variable << " float"
	 << "\nLOOKUP_TABLE default \n";
    for(int k = 1; k <= Nvtk_z ; ++k) {
	for(int j = 1; j <= Nvtk_y ; ++j) {
	    for(int i = 1; i <= Nvtk_x; ++i) {
		file << data(i,j) << " ";
		if ( !(i % 6) ) file << "\n";
	    }
	}
    }
        
    file.close();

}

void readFromFile_DX_3D(ScalarField3D& scf1, ScalarField3D& scf2, ScalarField3D& scf3, 
			int number, string filename, int w) 
{   
  //  int w = 5;
  ostringstream extnumber;
    extnumber.width(w);
    extnumber.fill('0');
    extnumber << number;
    string snumber = extnumber.str();
    filename += snumber;
    ifstream file (filename.c_str());

    if (!file) {
	cout << "\n !Can't open file \" " << filename << "\"" << endl;
    } else {
      //      cout << "\n Opening file \" " << filename << "\"" << endl;
    }
    
    char line[256];
/// 
/// by the time, just jump the first 8 lines
///	
    for(int ii = 0; ii < 9; ++ii) {
	file.getline(line,256);
	//	cout << " ii = " << ii << "\t" << line << endl;
    }

    file.precision(10);
    file.setf(ios_base::fixed,ios_base::floatfield);
///
/// Here, scf1 and scf2 are supposed to be of type ScalarField2D
///
    for(int k = 0; k <= scf1.ubound(thirdDim); ++k) 
      for(int i = 0; i <= scf1.ubound(firstDim); ++i) 
	for(int j = 0; j <= scf1.ubound(secondDim); ++j)
	  file >> scf1(i,j,k) >> scf2(i,j,k) >> scf3(i,j,k);
    file.close();
}


